#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <algorithm>

using namespace std;

// 函数类型枚举
enum FunctionType {
    LINEAR = 1,     // 一次函数
    QUADRATIC = 2,  // 二次函数
    CUBIC = 3       // 三次函数
};

// 拟合方法枚举
enum FitMethod {
    LEAST_SQUARES,  // 最小二乘法
    GRADIENT_DESCENT // 梯度下降法
};

// 数据点结构
struct DataPoint {
    double x;
    double y;
};

// 曲线拟合类
class CurveFitter {
private:
    vector<DataPoint> data;
    vector<double> fittedParams;
    double rSquared;
    FunctionType funcType;
    FitMethod fitMethod;

public:
    // 添加公共方法获取数据
    vector<DataPoint> getData() const {
        return data;
    }
    
    // 生成数据
    void generateData(FunctionType type, vector<double> params, int numPoints, double noiseLevel, const string& filename) {
        funcType = type;
        data.clear();
        
        // 初始化随机数生成器
        srand(time(NULL));
        
        // 生成数据点
        for (int i = 0; i < numPoints; ++i) {
            DataPoint point;
            // 生成-100到100之间的随机x值
            point.x = ((double)rand() / RAND_MAX) * 200.0 - 100.0;
            
            // 根据函数类型计算y值
            switch (type) {
                case LINEAR:
                    if (params.size() < 2) throw invalid_argument("线性函数需要2个参数");
                    point.y = params[0] * point.x + params[1];
                    break;
                case QUADRATIC:
                    if (params.size() < 3) throw invalid_argument("二次函数需要3个参数");
                    point.y = params[0] * pow(point.x, 2) + params[1] * point.x + params[2];
                    break;
                case CUBIC:
                    if (params.size() < 4) throw invalid_argument("三次函数需要4个参数");
                    point.y = params[0] * pow(point.x, 3) + params[1] * pow(point.x, 2) + params[2] * point.x + params[3];
                    break;
            }
            
            // 添加噪声
            double noise = ((double)rand() / RAND_MAX) * 2.0 - 1.0;
            point.y += noise * noiseLevel;
            data.push_back(point);
        }
        
        // 保存数据到文件
        saveData(filename);
    }
    
    // 从文件加载数据
    void loadData(const string& filename) {
        data.clear();
        
        ifstream file(filename);
        if (!file.is_open()) throw runtime_error("无法打开文件: " + filename);
        
        double x, y;
        while (file >> x >> y) {
            data.push_back({x, y});
        }
        
        file.close();
    }
    
    // 保存数据到文件
    void saveData(const string& filename) const {
        ofstream file(filename);
        if (!file.is_open()) throw runtime_error("无法创建文件: " + filename);
        
        for (vector<DataPoint>::const_iterator it = data.begin(); it != data.end(); ++it) {
            file << it->x << " " << it->y << endl;
        }
        
        file.close();
    }
    
    // 使用最小二乘法拟合
    void fitLeastSquares(FunctionType type) {
        funcType = type;
        fitMethod = LEAST_SQUARES;
        
        int n = data.size();
        if (n == 0) throw runtime_error("没有数据点可供拟合");
        
        switch (type) {
            case LINEAR:
                fitLinearLeastSquares();
                break;
            case QUADRATIC:
                fitQuadraticLeastSquares();
                break;
            case CUBIC:
                fitCubicLeastSquares();
                break;
        }
        
        calculateRSquared();
    }
    
    // 使用梯度下降法拟合
    void fitGradientDescent(FunctionType type, double learningRate = 0.0000001, int maxIterations = 100000, double tolerance = 1e-6) {
        funcType = type;
        fitMethod = GRADIENT_DESCENT;
        
        int n = data.size();
        if (n == 0) throw runtime_error("没有数据点可供拟合");
        
        // 初始化参数
        int paramCount;
        switch (type) {
            case LINEAR:
                paramCount = 2;
                break;
            case QUADRATIC:
                paramCount = 3;
                break;
            case CUBIC:
                paramCount = 4;
                break;
            default:
                throw invalid_argument("不支持的函数类型");
        }
        
        fittedParams.resize(paramCount);
        for (int i = 0; i < paramCount; ++i) {
            fittedParams[i] = 0.0; // 初始化为0
        }
        
        // 梯度下降迭代
        vector<double> oldParams(paramCount);
        for (int iter = 0; iter < maxIterations; ++iter) {
            // 保存旧参数
            for (int i = 0; i < paramCount; ++i) {
                oldParams[i] = fittedParams[i];
            }
            
            // 计算梯度
            vector<double> gradient = computeGradient();
            
            // 更新参数
            for (int i = 0; i < paramCount; ++i) {
                fittedParams[i] -= learningRate * gradient[i];
            }
            
            // 检查收敛
            double diff = 0.0;
            for (int i = 0; i < paramCount; ++i) {
                diff += pow(fittedParams[i] - oldParams[i], 2);
            }
            
            if (sqrt(diff) < tolerance) {
                break;
            }
        }
        
        calculateRSquared();
    }
    
    // 计算梯度
    vector<double> computeGradient() const {
        int n = data.size();
        int paramCount = fittedParams.size();
        vector<double> gradient(paramCount, 0.0);
        
        for (vector<DataPoint>::const_iterator it = data.begin(); it != data.end(); ++it) {
            double error = predict(it->x) - it->y;
            
            switch (funcType) {
                case LINEAR:
                    gradient[0] += error * it->x;
                    gradient[1] += error;
                    break;
                case QUADRATIC:
                    gradient[0] += error * pow(it->x, 2);
                    gradient[1] += error * it->x;
                    gradient[2] += error;
                    break;
                case CUBIC:
                    gradient[0] += error * pow(it->x, 3);
                    gradient[1] += error * pow(it->x, 2);
                    gradient[2] += error * it->x;
                    gradient[3] += error;
                    break;
            }
        }
        
        // 取平均
        for (int i = 0; i < paramCount; ++i) {
            gradient[i] /= n;
        }
        
        return gradient;
    }
    
    // 预测函数
    double predict(double x) const {
        switch (funcType) {
            case LINEAR:
                if (fittedParams.size() < 2) throw runtime_error("参数不足");
                return fittedParams[0] * x + fittedParams[1];
            case QUADRATIC:
                if (fittedParams.size() < 3) throw runtime_error("参数不足");
                return fittedParams[0] * pow(x, 2) + fittedParams[1] * x + fittedParams[2];
            case CUBIC:
                if (fittedParams.size() < 4) throw runtime_error("参数不足");
                return fittedParams[0] * pow(x, 3) + fittedParams[1] * pow(x, 2) + fittedParams[2] * x + fittedParams[3];
            default:
                throw invalid_argument("不支持的函数类型");
        }
    }
    
    // 计算R²值
    void calculateRSquared() {
        int n = data.size();
        if (n == 0) throw runtime_error("没有数据点可供计算R²");
        
        // 计算平均值
        double yMean = 0.0;
        for (vector<DataPoint>::const_iterator it = data.begin(); it != data.end(); ++it) {
            yMean += it->y;
        }
        yMean /= n;
        
        // 计算总平方和
        double ssTotal = 0.0;
        for (vector<DataPoint>::const_iterator it = data.begin(); it != data.end(); ++it) {
            ssTotal += pow(it->y - yMean, 2);
        }
        
        // 计算残差平方和
        double ssResidual = 0.0;
        for (vector<DataPoint>::const_iterator it = data.begin(); it != data.end(); ++it) {
            double yPred = predict(it->x);
            ssResidual += pow(it->y - yPred, 2);
        }
        
        // 计算R²
        rSquared = 1.0 - (ssResidual / ssTotal);
    }
    
    // 获取拟合参数
    vector<double> getFittedParams() const {
        return fittedParams;
    }
    
    // 获取R²值
    double getRSquared() const {
        return rSquared;
    }
    
    // ASCII可视化
    void visualizeASCII(int width = 80, int height = 20) const {
        if (data.empty()) {
            cout << "没有数据可显示" << endl;
            return;
        }
        
        // 找到数据范围
        double xMin = data[0].x, xMax = data[0].x;
        double yMin = data[0].y, yMax = data[0].y;
        
        for (vector<DataPoint>::const_iterator it = data.begin(); it != data.end(); ++it) {
            if (it->x < xMin) xMin = it->x;
            if (it->x > xMax) xMax = it->x;
            if (it->y < yMin) yMin = it->y;
            if (it->y > yMax) yMax = it->y;
        }
        
        // 扩展范围10%，使图像更美观
        double xRange = xMax - xMin;
        double yRange = yMax - yMin;
        xMin -= xRange * 0.1;
        xMax += xRange * 0.1;
        yMin -= yRange * 0.1;
        yMax += yRange * 0.1;
        
        // 创建画布
        vector<string> canvas(height, string(width, ' '));
        
        // 绘制坐标轴
        int xAxisY = static_cast<int>((yMax - 0) / (yMax - yMin) * (height - 1));
        int yAxisX = static_cast<int>((0 - xMin) / (xMax - xMin) * (width - 1));
        
        if (xAxisY >= 0 && xAxisY < height) {
            for (int x = 0; x < width; ++x) {
                canvas[xAxisY][x] = '-';
            }
        }
        
        if (yAxisX >= 0 && yAxisX < width) {
            for (int y = 0; y < height; ++y) {
                canvas[y][yAxisX] = '|';
            }
        }
        
        // 绘制原点
        if (xAxisY >= 0 && xAxisY < height && yAxisX >= 0 && yAxisX < width) {
            canvas[xAxisY][yAxisX] = '+';
        }
        
        // 绘制数据点
        for (vector<DataPoint>::const_iterator it = data.begin(); it != data.end(); ++it) {
            int xPixel = static_cast<int>((it->x - xMin) / (xMax - xMin) * (width - 1));
            int yPixel = static_cast<int>((yMax - it->y) / (yMax - yMin) * (height - 1));
            
            if (xPixel >= 0 && xPixel < width && yPixel >= 0 && yPixel < height) {
                canvas[yPixel][xPixel] = '*';
            }
        }
        
        // 绘制拟合曲线
        for (int xPixel = 0; xPixel < width; ++xPixel) {
            double x = xMin + (xMax - xMin) * xPixel / (width - 1);
            double y = predict(x);
            int yPixel = static_cast<int>((yMax - y) / (yMax - yMin) * (height - 1));
            
            if (yPixel >= 0 && yPixel < height) {
                if (canvas[yPixel][xPixel] == ' ') {
                    canvas[yPixel][xPixel] = '+';
                } else if (canvas[yPixel][xPixel] == '*') {
                    canvas[yPixel][xPixel] = 'x';
                }
            }
        }
        
        // 输出画布
        cout << endl;
        for (vector<string>::const_iterator it = canvas.begin(); it != canvas.end(); ++it) {
            cout << *it << endl;
        }
        
        // 输出拟合信息
        cout << endl << "拟合参数: ";
        for (size_t i = 0; i < fittedParams.size(); ++i) {
            cout << static_cast<char>('a' + i) << " = " << fittedParams[i] << "  ";
        }
        cout << endl << "R² = " << rSquared << endl;
    }
    
private:
    // 一次函数最小二乘法拟合
    void fitLinearLeastSquares() {
        int n = data.size();
        double sumX = 0.0, sumY = 0.0, sumXY = 0.0, sumX2 = 0.0;
        
        for (vector<DataPoint>::const_iterator it = data.begin(); it != data.end(); ++it) {
            sumX += it->x;
            sumY += it->y;
            sumXY += it->x * it->y;
            sumX2 += it->x * it->x;
        }
        
        double denominator = n * sumX2 - sumX * sumX;
        if (denominator == 0) {
            throw runtime_error("无法拟合线性函数: 分母为零");
        }
        
        double a = (n * sumXY - sumX * sumY) / denominator;
        double b = (sumY * sumX2 - sumX * sumXY) / denominator;
        
        fittedParams.resize(2);
        fittedParams[0] = a;
        fittedParams[1] = b;
    }
    
    // 二次函数最小二乘法拟合
    void fitQuadraticLeastSquares() {
        int n = data.size();
        double sumX = 0.0, sumX2 = 0.0, sumX3 = 0.0, sumX4 = 0.0;
        double sumY = 0.0, sumXY = 0.0, sumX2Y = 0.0;
        
        for (vector<DataPoint>::const_iterator it = data.begin(); it != data.end(); ++it) {
            double x = it->x;
            double x2 = x * x;
            double x3 = x2 * x;
            double x4 = x3 * x;
            
            sumX += x;
            sumX2 += x2;
            sumX3 += x3;
            sumX4 += x4;
            sumY += it->y;
            sumXY += x * it->y;
            sumX2Y += x2 * it->y;
        }
        
        // 构建并求解正规方程组
        double det = n*sumX2*sumX4 + 2*sumX*sumX2*sumX3 - sumX2*sumX2*sumX2 - n*sumX3*sumX3 - sumX*sumX*sumX4;
        if (det == 0) {
            throw runtime_error("无法拟合二次函数: 行列式为零");
        }
        
        double a = (sumY*sumX2*sumX4 + sumXY*sumX3*sumX2 + sumX2Y*sumX*sumX3 - sumY*sumX3*sumX3 - sumXY*sumX*sumX4 - sumX2Y*sumX2*sumX2) / det;
        double b = (n*sumXY*sumX4 + sumY*sumX3*sumX2 + sumX2Y*sumX2*sumX - n*sumX2Y*sumX3 - sumXY*sumX2*sumX2 - sumY*sumX*sumX4) / det;
        double c = (n*sumX2*sumX2Y + sumX*sumX2*sumXY + sumX*sumX3*sumY - n*sumX3*sumXY - sumX*sumX*sumX2Y - sumX2*sumX2*sumY) / det;
        
        fittedParams.resize(3);
        fittedParams[0] = a;
        fittedParams[1] = b;
        fittedParams[2] = c;
    }
    
    // 三次函数最小二乘法拟合
    void fitCubicLeastSquares() {
        int n = data.size();
        double sumX = 0.0, sumX2 = 0.0, sumX3 = 0.0, sumX4 = 0.0, sumX5 = 0.0, sumX6 = 0.0;
        double sumY = 0.0, sumXY = 0.0, sumX2Y = 0.0, sumX3Y = 0.0;
        
        for (vector<DataPoint>::const_iterator it = data.begin(); it != data.end(); ++it) {
            double x = it->x;
            double x2 = x * x;
            double x3 = x2 * x;
            double x4 = x3 * x;
            double x5 = x4 * x;
            double x6 = x5 * x;
            
            sumX += x;
            sumX2 += x2;
            sumX3 += x3;
            sumX4 += x4;
            sumX5 += x5;
            sumX6 += x6;
            
            sumY += it->y;
            sumXY += x * it->y;
            sumX2Y += x2 * it->y;
            sumX3Y += x3 * it->y;
        }
        
        // 构建正规方程组的系数矩阵和常数向量
        // 修改：使用正确的变量名sumX, sumX2, sumX3等
        vector<vector<double>> A = {
            {static_cast<double>(n), sumX, sumX2, sumX3},
            {sumX, sumX2, sumX3, sumX4},
            {sumX2, sumX3, sumX4, sumX5},
            {sumX3, sumX4, sumX5, sumX6}
        };
        
        // 修改：使用正确的变量名sumY, sumXY, sumX2Y, sumX3Y
        vector<double> b = {sumY, sumXY, sumX2Y, sumX3Y};
        
        // 使用高斯消元法求解方程组
        vector<double> params = gaussianElimination(A, b);
        
        fittedParams = params;
    }
    
    // 高斯消元法求解线性方程组
    vector<double> gaussianElimination(vector<vector<double>> A, vector<double> b) {
        int n = A.size();
        
        // 前向消元
        for (int i = 0; i < n; ++i) {
            // 寻找主元
            int max_row = i;
            for (int j = i + 1; j < n; ++j) {
                if (fabs(A[j][i]) > fabs(A[max_row][i])) {
                    max_row = j;
                }
            }
            
            // 交换行
            swap(A[i], A[max_row]);
            swap(b[i], b[max_row]);
            
            // 检查主元是否为零
            if (fabs(A[i][i]) < 1e-10) {
                throw runtime_error("矩阵奇异，无法求解");
            }
            
            // 归一化当前行
            double pivot = A[i][i];
            for (int j = i; j < n; ++j) {
                A[i][j] /= pivot;
            }
            b[i] /= pivot;
            
            // 消元
            for (int j = i + 1; j < n; ++j) {
                double factor = A[j][i];
                for (int k = i; k < n; ++k) {
                    A[j][k] -= factor * A[i][k];
                }
                b[j] -= factor * b[i];
            }
        }
        
        // 回代
        vector<double> x(n);
        for (int i = n - 1; i >= 0; --i) {
            x[i] = b[i];
            for (int j = i + 1; j < n; ++j) {
                x[i] -= A[i][j] * x[j];
            }
        }
        
        return x;
    }
};

// 命令行参数解析
struct CommandLineArgs {
    string command;         // "gendat" 或 "fitting"
    string filename;        // 数据文件名
    FunctionType funcType;  // 函数类型
    vector<double> params;  // 函数参数
    FitMethod fitMethod;    // 拟合方法
    bool visualize;         // 是否可视化
};

// 解析命令行参数
CommandLineArgs parseArgs(int argc, char* argv[]) {
    CommandLineArgs args;
    
    if (argc < 2) {
        throw invalid_argument("请指定命令 (gendat 或 fitting)");
    }
    
    args.command = argv[1];
    
    if (args.command == "gendat") {
        if (argc < 5) {
            throw invalid_argument("gendat 命令需要至少4个参数: gendat filename t a b [c] [d]");
        }
        
        args.filename = argv[2];
        args.funcType = static_cast<FunctionType>(stoi(argv[3]));
        
        if (args.funcType < LINEAR || args.funcType > CUBIC) {
            throw invalid_argument("函数类型必须是1(线性)、2(二次)或3(三次)");
        }
        
        int requiredParams = (args.funcType == LINEAR) ? 2 : (args.funcType == QUADRATIC) ? 3 : 4;
        if (argc < 4 + requiredParams) {
            throw invalid_argument("函数类型 " + to_string(args.funcType) + " 需要 " + to_string(requiredParams) + " 个参数");
        }
        
        for (int i = 0; i < requiredParams; ++i) {
            args.params.push_back(stod(argv[4 + i]));
        }
        
        args.visualize = false;
    }
    else if (args.command == "fitting") {
        if (argc < 5) {
            throw invalid_argument("fitting 命令需要4个参数: fitting filename m t");
        }
        
        args.filename = argv[2];
        
        string method = argv[3];
        if (method == "ls") {
            args.fitMethod = LEAST_SQUARES;
        } else if (method == "gd") {
            args.fitMethod = GRADIENT_DESCENT;
        } else {
            throw invalid_argument("拟合方法必须是 'ls' (最小二乘法) 或 'gd' (梯度下降法)");
        }
        
        args.funcType = static_cast<FunctionType>(stoi(argv[4]));
        if (args.funcType < LINEAR || args.funcType > CUBIC) {
            throw invalid_argument("函数类型必须是1(线性)、2(二次)或3(三次)");
        }
        
        args.visualize = (argc > 5 && string(argv[5]) == "--visualize");
    }
    else {
        throw invalid_argument("未知命令: " + args.command + "，支持的命令: gendat, fitting");
    }
    
    return args;
}

int main(int argc, char* argv[]) {
    try {
        CommandLineArgs args = parseArgs(argc, argv);
        CurveFitter fitter;
        
        if (args.command == "gendat") {
            // 生成数据
            int numPoints = 100;  // 默认100个点
            double noiseLevel = 0.1;  // 默认噪声水平
            
            // 可选参数处理
            for (int i = 5 + args.params.size(); i < argc; i += 2) {
                string option = argv[i];
                if (option == "--num_points" && i + 1 < argc) {
                    numPoints = stoi(argv[i + 1]);
                } else if (option == "--noise_level" && i + 1 < argc) {
                    noiseLevel = stod(argv[i + 1]);
                } else {
                    throw invalid_argument("未知选项: " + option);
                }
            }
            
            fitter.generateData(args.funcType, args.params, numPoints, noiseLevel, args.filename);
            cout << "数据已生成并保存到 " << args.filename << endl;
        }
        else if (args.command == "fitting") {
            // 加载数据
            fitter.loadData(args.filename);
            // 使用公共接口获取数据大小
            cout << "已加载 " << args.filename << "，包含 " << fitter.getData().size() << " 个数据点" << endl;
            
            // 执行拟合
            if (args.fitMethod == LEAST_SQUARES) {
                fitter.fitLeastSquares(args.funcType);
            } else {
                fitter.fitGradientDescent(args.funcType);
            }
            
            // 输出结果
            vector<double> params = fitter.getFittedParams();
            cout << "拟合结果:" << endl;
            for (size_t i = 0; i < params.size(); ++i) {
                cout << static_cast<char>('a' + i) << " = " << params[i] << endl;
            }
            cout << "R² = " << fitter.getRSquared() << endl;
            
            // 可视化
            if (args.visualize) {
                fitter.visualizeASCII();
            }
        }
    }
    catch (const exception& e) {
        cerr << "错误: " << e.what() << endl;
        return 1;
    }
    return 0;
    }